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A NONLINEAR  FILTERING  TECHNIQUE  FOR  DIGITIZED 
IMAGES  DEGRADED  BY  FILM-GRAIN  NOISE 

I.  INTRODUCTION  AND  NOTATION 


One  kind  of  image  degradation  that  often  occurs  is  film-grain  noise,  which  is  produced  by 
a photographic  emulsion  during  the  process  of  recording  an  image  on  film.  An  extensive 
account  of  the  mechanisms  and  behavior  of  this  type  of  noise  can  be  found  in  Naderi  (ll, 
James  [2],  and  included  references.  This  report  investigates  the  use  of  digital  image-processing 
techniques  to  suppress  this  grain  noise  by  implementing  a particular  class  of  nonlinear  filters. 
The  emphasis  here  is  on  processing  an  image  in  terms  of  the  "optical  density"  variable,  which  is 
a logarithmic  (point-by-point)  function  of  the  light  intensity  transmitted  by  a film  negative, 
because  the  grain  noise  is  of  the  convenient  additive  form  in  this  variable. 

One  troublesome  feature  of  film-grain  noise  is  the  fact  that  it  is  signal  dependent,  i.e.,  its 
root-mean-square  (rms)  value  varies  with  the  local  average  optical  density  of  the  film.  The  con- 
ventional "optimal"  Wiener  filter  for  suppressing  this  noise  ignores  its  signal  dependence  ID. 
The  type  of  nonlinear  filter  applied  here,  however,  takes  the  first-order  effects  of  this  signal 
dependence  into  account,  and  experimentally  shows  some  improvement  over  the  conventional 
Wiener  filter.  This  nonlinear  filter  is  based  on  Bayesian  probability  concepts,  and  it  results 
from  a perturbation  analysis  of  a linear  filtering  problem.  Hence,  considerable  background  of 
relevant  linear  filtering  is  also  presented. 


Authorities  disagree  over  the  relative  importance  of  grain-noise  signal  dependence  and  the 
difficulties  inherent  in  the  nonlinear  relation  between  image  intensity  and  film  density  13,4).  It 
happens  that  the  class  of  nonlinear  filters  developed  here  can  also  be  used  to  account  for  such 
nonlinearities  in  images  of  sufficiently  low  contrast,  so  this  application  is  also  explored  to  some 
extent.  This  particular  nonlinearity  turns  out  to  have  rather  severe  effects  in  practice,  however, 
and  the  validity  of  the  perturbation  analysis  used  here  is  somewhat  dubious  in  this  context. 


Unless  otherwise  indicated,  lower-case  letters  in  the  following  material  denote  (real) 
column  vectors  or  scalars.  Matrices  are  denoted  by  capital  Roman  letters.  A ^denotes  the  tran- 
spose of  a matrix  A,  and  tr(A)  its  trace.  It  will  also  be  convenient  to  manipulate  three-way 
matrices,  which  are  denoted  by  capital  Greek  letters.  For  continuity  of  notation,  the  following 
definitions  are  adopted  for  such  a matrix  P,  with  vector  x and  matrices  A and  B of  compatible 
dimensions,  and  with  repeated  indices  denoting  summation: 


(I  x) , y “ r, 

^AP)  , ” A,„  P„  y* 

D B),  y*  P,  y„B„* 
(AxO,y»  - A,p(^ 


(matrix) 

(three-way  matrix) 
(three-way  matrix) 
(three-way  matrix) 
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(three-way  matrix) 

“ r*  y, 

(three-way  matrix) 

(Tr(r)iy-r„,„ 

(column  vector) 

ITr(rr)i,y-  rx,.,r„y. 

(matrix). 

Also,  r is  called  symmetric  if  r — T'  — T"  - F 7- 


II.  BACKGROUND 
A.  Linear  Filtering 

Because  of  its  virtually  complete  lack  of  spatial  correlation,  grain  noise  can  always  be 
reduced  by  local  spatial  averaging,  provided  of  course  that  the  averaging  is  done  over  an  area 
whose  dimensions  are  large  compared  to  that  of  a typical  grain  (~1  nm).  Perhaps  the  simplest 
procedure  of  this  sort  is  to  assign  the  optical  density 

— ^ zir.9)dedr 

Jo  Jo 

to  each  point  of  the  processed  image,  where  (r.O)  denotes  a polar  coordinate  system  centered  at 
the  point  of  interest,  z is  the  recorded  optical  density,  and  R is  some  small  radius.  In  other 
words,  the  density  at  each  point  in  the  processed  image  is  merely  the  average  of  the  observed 
density  in  a small  circle  of  radius  R centered  at  that  point.  In  digital  image  processing,  of 
course,  the  integrals  are  approximated  by  discrete  sums. 

The  averaging  nature  of  this  procedure  means  that  it  achieves  grain  noise  suppression 
only  at  the  cost  of  a loss  of  resolution.  Some  penalty  of  this  sort  is  always  unavoidable,  but 
more  sophisticated  'filtering''  methods  can  achieve  greater  grain-noise  suppression  for  a given 
reduction  in  resolution.  A common  type  of  procedure  is  linear  filtering,  which  corresponds  to 
weighted  spatial  averaging.  In  its  most  general  form,  this  means  (in  rectangular  coordinates) 
that  the  density  6 at  a point  (x,  y)  is  computed  as 

b(x,y)  - ;;  Vi{x,  y,u,\)z(u,\)  du  d\, 

where  the  region  of  integration  is  the  entire  image.  In  practice,  the  weighting  function 
W(x,  y,u,\)  is  almost  always  chosen  so  that  its  values  depend  only  on  x - uand  y - v (except 
near  the  e'tge  of  the  image).  The  computed  density  in  this  case  can  be  interpreted  as  the  two- 
dimensir  .j  convolution  of  the  weighting  function  with  the  observed  image;  that  is, 

<»(x,  y)  - ;;  W(0,0.u -x,v-y)  z(u,\)du  d\. 

This  can  be  expressed  more  concisely  as 

h(x,  y)  z(x.  y)  dx  dy, 

where  x and  y are  now  local  coordinates  centered  at  the  point  whose  filtered  density  is  being 
computed,  and  where  h(x,  y)  - W(0,0,x,  y). 

Clearly,  the  local  spatial  averaging  method  described  above  is  a particular  kind  of  linear 
filter,  with 
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It  is  generally  regarded  as  reasonable  to  approximate  grain  noise,  within  a uniformly 
exposed  region  at  least,  as  spatially  uniform  additive  white  noise  because  of  the  absence  of  spa- 
tial correlation.  Without  going  into  details,  this  means  that  the  relative  grain-noise  suppression 
for  a weighting  function  h (x,  y)  can  be  expressed  as 


where 

k ••  h(x,y)dxdy  (a  normalization  constant) 

a,  — rms  noise-induced  density  fluctuation  in  processed  image  (using  h) 
o-g  ~ corresponding  rms  fluctuation  with  a uniform  window  of  unit  area. 


In  digital  image  processing,  the  pixel  size  is  generally  used  as  the  unit  area.  This  result  enables 
us  to  And,  for  example,  the  relative  sizes  of  R,  cr,  and  a in  the  preceding  three  filters  that  yield 
equal  grain-noise  suppression,  namely 


The  form  of  these  three  parameters  is  such  that  increasing  them  increases  the  degree  of  grain- 
noise  suppression,  and  decreases,  at  the  same  time,  the  resolution.  The  differing  tradeoffs 
between  resolution  and  noise  suppression  can  now  be  compared  by  plotting  the  responses  of 
these  three  filters  to  a semi-infinite  step  for  equal  grain-noise  suppression.  The  results  are  sum- 
marized in  Fig.  2,  and  are  independent  of  the  step  size  and  degree  of  noise  suppression  except  ^ 

for  scale  changes.  The  differences  are  not  dramatic,  but  this  example  does  illustrate  how  this 
tradeoff  can  be  affected  by  the  type  of  filter  used.  For  instance,  if  one  adopts  as  a criterion  of 

resolution  the  distance  required  for  the  response  to  rise  from  one-quarter  to  three-quarters  of  t 

its  maximum  value,  then  for  the  orientation  shown,  the  exponential  window  is  about  20% 
better  than  the  uniform  circular  window  with  the  same  degree  of  grain-noise  suppression.  It  is 
even  better  for  other  orientations. 

B.  Bayesian  Filtering 

Another  approach  to  this  problem  is  to  postulate  or  otherwise  determine  an  a priori  proba- 
bility distribution  for  the  set  of  possible  original  images  and  then  compute  its  conditional  proba-  j 

bility  distribution,  given  the  observed  noise-corrupted  image.  In  terms  of  a spatially  discrete 
image,  whose  pixels  are  regarded  for  convenience  as  continuously  variable  in  optical  density, 
this  typically  means  specifying  the  probability  densities  p(x)  and  p(zjx,),  with  i — I,  ....  N-, 
where  x is  an  N-dimensional  real  vector  whose  components  correspond  to  the  optical  densitie.s 

of  the  individual  pixels,  and  where  z,  is  the  observed  optical  density  of  the  / th  pixel.  The  i 

desired  posterior  probability  density  p(x/Z|,  ■ • ■ , zs)  can  then  be  determined  in  principle  by  | 

using  the  Bayes  rule.  The  number  N of  pixels  in  an  image  is  typically  in  the  hundreds  of  ! 

thousands,  however,  and  such  a procedure  is  not  computationally  feasible  at  present  except  in 
special  cases.  Another  drawback  to  this  approach  is  that,  whereas  p(zjx)  is  often  definitely 

known,  the  prior  probability  density  p(x)  usually  is  not,  and  is  fairly  arbitrary.  I 
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Rig.  2 — Responses  to  semi-inlinile  step  for  three  two-dimensional  filters  J 

with  equal  grain-noise  suppression 


With  the  Bayesian  approach,  one  is  usually  content  with  finding  either  the  mean  or  mode 
of  the  posterior  probability  density  as  a good  estimate  of  the  original  image,  given  the  observed 
image.  An  important  property  in  this  regard  is  that  if 

p(zjx)  - X,-\-  VD,  (1) 

and  if  p{x)  and  p(w^,  • ■ • . w^v)  are  both  multivariate  Gaussian  probability  densities  with  a 

priori  known  parameters  (mean  and  covariance),  then  the  posterior  probability  density  p(x/z)  I 

is  also  multivariate  Gaussian.  This  means  that  its  mean  and  mode  are  identical.  Furthermore, 
the  value  of  this  conditional  mean  for  each  pixel  is  a linear  combination  of  the  observed  optical 
densities  of  the  pixels.  Hence,  the  process  of  computing  the  conditional  mean,  which  we  will 
denote  as  Bayesian  Altering  here,  falls  into  the  category  of  linear  Alters  described  in  the  preced- 
ing section  under  these  particular  conditions  of  additive  Gaussian  measurement  noise  and  a 
Gaussian  prior.  Incidentally,  the  Bayes  Alter  in  this  context  is  also  the  same  as  the  Wiener 
Alter,  which  is  optimal  in  the  least  squares  sense. 

Because  of  this  property,  one  way  of  constructing  linear  Alters  for  grain-noise  suppression 
is  to  treat  the  noise  as  additive  and  Gaussian  and  to  use  the  corresponding  Bayes  Alter  for  some 
Gaussian  prior  distribution  that  is  roughly  consistent  with  the  situation  at  hand.  The  rationale 
for  this  construction  is  that  the  resulting  linear  Alter  is  optimal  in  at  least  some  sense  for  a simi- 
lar situation,  and  therefore  hopefully  good  for  the  image  actually  being  processed.  However, 
this  procedure  requires  that  the  grain-noise  magnitude  be  approximated  as  independent  of  the 
optical  density  of  the  underlying  image,  which  is  considerably  at  variance  with  experimental  evi- 
dence (5).  Much  of  the  emphasis  in  this  report  is  on  modifying  the  Bayes  Alter  by  perturbation 
methods  to  account  more  accurately  for  dependence,  of  the  grain  noise  on  optical  density.  The 
resulting  Alters  are  no  longer  linear,  however,  and  cannot  be  interpreted  as  weighted  spatial 
averaging. 
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As  a concession  to  computational  tractability,  this  Bayesian  filtering  approach  can  be 
applied  to  the  pixel  rows  and  columns  of  a discrete  image  individually,  resulting  in  a one- 
dimensional filter  for  each  row  or  column  (1).  Each  such  filter  is  based  on  an  a priori  joint  pro- 
bability distribution  for  the  optical  density  values  of  only  the  pixels  in  the  row  or  column  in 
question.  This  prior  distribution  is  chosen  to  be  consistent  with  the  main  features  of  the  image, 
such  as  resolution  and  rms  optical  density  fluctuation.  The  filtering  then  amounts  to  computing 
the  posterior  mean  of  the  pixel  densities  in  that  row  or  column  only,  given  the  noise-corrupted 
observations  thereof.  The  grain  noise,  being  statistically  independent  for  each  pixel,  is  the 
same  as  in  the  two-dimensional  analysis.  Because  it  treats  each  row  and  column  individually, 
this  one-dimensional  analysis  does  not  address  the  possibility  of  correlations  between  adjacent 
rows  or  columns  in  the  original  image,  which  typically  exist  to  a high  degree,  nor  the  fact  that 
rows  and  columns  intersect.  Realistically,  then,  the  use  of  such  methods  in  two-dimensional 
image  processing  can  be  only  partly  justified  on  Bayesian  grounds.  Also,  the  further 
simplification  is  often  made  of  using  the  same  prior  distribution,  and  hence  the  same  filter,  for 
each  row  and  column.  This  is  a reasonable  approximation  in  view  of  the  other  inaccuracies 
already  incurred. 

If  the  entire  image  is  processed,  these  one-dimensional  filters  are  then  combined  by  some 
other  rationale.  One  procedure,  described  by  Naderi  (1),  is  to  filter  each  row  and  column  of  the 
observed  image  separately,  and  then  assign  to  each  pixel  in  the  processed  image  the  optical  den- 
sity that  is  the  average  of  the  filtered  values  for  that  pixel  in  the  row  and  column  to  which  it 
belongs.  Another  method  is  to  filter  each  row  separately  and  then  filter  the  columns  of  the 
resulting  image  individually.  If  the  row  and  column  filters  are  linear  and  if  their  results  are  com- 
bined linearly  in  processing  the  two-dimensional  image,  then  the  composite  filtering  operation 
is  a linear  filter  of  the  type  described  in  the  preceding  section,  and  the  weighting  function  of  the 
corresponding  weighted  spatial  average  can  be  examined  directly  as  a guide  in  combining  the 
constituent  one-dimensional  filters.  For  linear  row  and  column  filters,  both  procedures  just 
mentioned  are  linear  combinations  of  this  sort.  A special  case  of  particular  interest  occurs 
when  the  row  and  column  filters  all  correspond  to  a single  one-dimensional  weighting  function 
of  the  spatially  invariant  Gaussian  form 

.jL 

h(x)  - k e k,  (T  are  known  constants 

except  near  the  edges.  This  case  is  noteworthy  because  if  these  one-dimensional  filters  are 
combined  in  the  second  way  described  above  (i.e.,  by  cascading  the  row  and  column  filtering 
operations),  then,  except  for  edge  effects,  the  resulting  two-dimensional  linear  filter  has  the 
"Gaussian  window"  type  of  weighting  function  described  in  the  preceding  section.  This  weight- 
ing function  is  desirable  because  of  its  isotropic  properties. 

Most  of  this  report  is  further  restricted  to  filters  composed  of  one-dimensional  Bayes 
filters  for  rows  and  columns  of  image  pixels.  Although  these  constituent  filters  are  generally 
nonlinear,  the  types  we  examine  are  derived  from  a perturbation  analysis  of  linear  filters,  so  the 
weighted  spatial  averaging  interpretation  of  the  corresponding  unperturbed  linear  filters  is  used 
as  a guide  in  combining  the  nonlinear  filters  into  two-dimensional  filters. 


I 


NRL  REPORT  8225 


I 


i 

I 


I 


C.  Recursive  Implementations 

It  is  sometimes  possible  to  reduce  the  computational  burden  of  Bayesian  filtering  by 
describing  the  optical  density  fluctuations  of  the  original  image  as  a Markov  process  and  exploit- 
ing this  Markov  property  to  construct  a recursive,  or  sequential,  implementation  of  the 
corresponding  Bayesian  filter.  This  approach  is  emphasized  in  this  report,  and  is  especially 
suited  to  the  type  of  one-dimensional  line  filters  Just  described. 

As  an  example,  consider  the  case  of  the  density-independent  grain  noise  of  Eq.  (1)  with  a 
Gaussian  prior  distribution  for  the  original  image  fluctuations  having  the  autocorrelation  func- 
tion 

+ (2) 

where  u,  , denotes  the  optical  density  of  the  i,  / th  pixel.  For  simplicity,  regard  all  random  vari- 
ables as  having  zero  mean,  which  introduces  no  loss  of  generality.  According  to  the  theory 
described  in  the  preceding  section,  one  could  aways  compute  the  (Gaussian)  conditional  proba- 
bility distribution  of  the  original  image,  given  the  observed  image  in  a single  "step"  by  defining 
a vector  x composed  of  all  the  image  pixels.  The  conditional  mean  x and  covariance  matrix  P 
of  this  vector  (which  completely  determine  this  posterior  distribution)  would  then  be  given  by 
the  well-known  equations  [6] 


M - -5-  M(I  -1-  M)-'M 
r 

(3) 

— Pz,  (filtered  image  vector) 

(4) 

where 

z " vector  of  observed  pixel  values 

M — covariance  matrix  of  a priori  distribution  of  x 

r * grain  noise  variance  (a  scalar  constant). 

If  the  image  is  an  N x N array  of  pixels,  however,  x,  x,  and  z here  are  vectors,  and  P and 
M are  x matrices.  For  realistic  values  of  N (often  512),  it  is  impossible  even  to  store  M 
or  P in  current  computers. 

According  to  Habibi  (7),  the  image  statistics  of  Eq.  (2)  can  be  obtained  as  the  output  of 
the  two-dimensional  Markov  process  with 

U,  j - p(u,-i.  , -I-  M,  y_,)  - p^M,_|  -I-  w,  !.  (5) 

where  [w,  /.  /,  y-0,  • • • , Nj  is  a set  of  independent  zero-mean  Gaussian  random  variables, 
each  with  variance  o-^(l  - p^,  in  which 

p - (6) 

It  is  apparent  from  the  form  of  Eq.  (5)  that  the  sequence  of  diagonals  in  the  pixel  array  can  be 
realized  as  a one-dimensional-vector  Markov  process  by  constructing  a state  vector  consisting  of 
pairs  of  diagonals.  For  convenience,  we  assume  a square  N x N image  and  make  all  the  state 
vectors  in  the  sequence  of  the  same  dimension  by  adding  fictitious  components  to  the  diagonals 
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beyond  the  edge  of  the  actual  image.  If  the  / th  augmented  diagonal  is  denoted  by  y,  (an  N 
vector),  then  this  is  equivalent  to  a Markov  process  with  the  linear  "dynamics"  of 


_ 0_ 

I__ 

Jl. 

_p  _ 

-pV 

P^/+2 

-(- 

“,  + 2 

/-I,  • • • , N-2. 


(7) 


where  the  w,  are  a sequence  of  independent  vector  zero-mean  Gaussian  random  variables,  each 
with  a diagonal  covariance  Q having  all  diagonal  elements  equal  to  <r^(l  - p^.  The  quantity  K, 
is  a matrix  of  the  form 


y 

0 10 0 T 

1 0^  i I 
0 \ 0 ' 

1 oil 

0 0 1 o-t- 

h — ' — " 


0 

0 


where  the  nonzero  diagonal  "block"  shown  in  detail  is  in  the  center  of  the  matrix.  The 
observed  pixel  densities  can  also  be  expressed  in  terms  of  this  sequence  of  state  vectors  as 


C,  - [OlH,] 


(8) 


where  i,  is  now  an  N-vector  of  observations  and  the  n,  are  independent,  identically  distributed 
zero-mean  vector  Gaussian  random  variables  with  diagonal  covariance  matrix  R,  all  of  whose 
diagonal  terms  are  equal  to  r.  The  matrices  H,  are  of  the  diagonal  form 


H, 


''^0  A 
'nV 


'o 


where  the  / nonzero  terms  are  in  the  center  of  the  matrix.  With  this  construction,  the  standard 
theory  of  Kalman  filtering  and  "linear  smoothing"  can  be  applied  (6)  to  the  realization  of  Eqs. 
(7)  and  (8)  to  obtain  the  posterior  probability  distribution  of  the  optical  densities  of  the  image 
pixels  according  to  the  following  formulas: 


F,x,  + ^ P,+|H4,({,+2  - H,+,F,x,) 

(9) 

M,  - M,H,^(H,M,H,''-t- 

(10) 
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M 


1 + 1 


B, 


F,P,F/+  Q 

X,  + P,F,^Mr+i(x,  + i - 

P,  - P,F,^Mr+',(M,^,  - 


F,x,); 

- 


(11) 

Xn  - x^ 

(12) 

/+'iF/Pi;  B/v“P^. 

(13) 

where 


h,-(o;h,^,],  Q 


_0_ 

1 

-LP. 

o’ 

1 Q 

and  F 


0 

-p>" 

This  equation  system  has  the  form  of  a forward  recursion,  consisting  of  Eqs. (9) -(1 1), followed 
by  a reverse  recursion  (Eqs.  (12)  and  (13)),  which  uses  the  results  of  the  forward  recursion. 
The  initial  values  Xj  and  Mj  are  specified  a priori.  The  interpretation  here  is  that  the  posterior 


distribution  of  the  composite  vector 


is  Gaussian,  with  mean  x,  and  covariance  mairi*  3,. 


y,+i 


Although  the  description  of  this  recursive  implementation  is  more  complicated,  the 
dimension  of  the  vectors  and  matrices  is  2N  instead  of  N^.  This  reduces  the  actual  computa- 
tional burden  enormously,  typically  by  a factor  of  100,000,  and  the  storage  requirements  to  an 
even  greater  extent.  However,  the  above  recursions  do  not  determine  the  "cross-correlations" 
Efy.y,)  for  |(— y|  >1.  These  features  of  the  posterior  distribution  are  not  usually  of  interest, 
and  in  any  event  can  be  computed  from  other  known  formulas  [8].  The  computational  require- 
ments of  this  recursive  procedure  are  still  formidable,  though  probably  within  the  realm  of  pos- 
sibility. Hence,  it  will  not  be  pursued  further  in  this  report. 

Another  particular  recursive  implementation,  which  will  be  pursued  in  greater  detail  later, 
is  for  a class  of  one-dimensional  Bayesian  line  (Iters.  The  one-dimensional  version  of  the  auto- 
correlation function  (2)  for  the  original  image  is 

E(u,h.„«,)  - (14) 

This  can  be  obtained,  for  the  zero-mean  Gaussian  case,  as  the  output  of  the  scalar  Markov 
process 

M,  I = /u,  -t-w, — 1,-  -,N  — 1 (15) 

in  which  the  w,  are  independent,  identically  distributed,  scalar,  Gaussian  random  variables  with 
variance  q,  and  where 

/-e-*  (16) 

and 

q — (r^{\  — e ^*).  (17) 

With  the  same  grain-noise  model  as  before  for  the  observed  image,  the  observed  optical  density 
z,  of  the  / th  pixel  in  such  a line  is  clearly 

z,-  u,  + v,;  1-0.  ■ • • . N (18) 

where  the  v,  are  random  variables  with  variance  r.  The  posterior  mean  and  variance  of  the  opt- 
ical densities  of  the  individual  pixels  are  then  given  by  the  scalar  version  of  the  recursions 
(9)  -(13),  which  is 
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+ , - fu,  + ^ (2,  + , - /m,);  m„  - 2„ 

m,r 
m,  + r 


y Kalman  filter 


"»,  + i - f^P,  + <7 


and 


fP, 


b,^  P,- 


"»,+i 

fP, 


"»,+i 


(w,+i  - /«,);  M/v  - 


(rW(  + i ~ ^/+i);  bf^  “ 


> linear  smoother 


(19) 

(20) 
(21) 

(22) 

(23) 


To  reiterate,  the  conditional  probability  density  of  «„  given  the  observed  image  (which  is  only  a 
line  here),  is  Gaussian  with  mean  u,  and  variance  b,. 


Except  near  the  edges,  the  p,  m and  b variables  are  essentially  constant  as  functions  of  /. 
This  permits  us  to  reduce  the  Bayes  filter  of  Eqs.  (19)  and  (22)  to  a one-dimensional  convolu- 
tion filter  of  the  form 


V 

I 


/-O 


hU  -j)z(j). 


(24) 


Working  out  the  details,  we  note  that  h in  Eq.  (24)  is  given  by 


h(n) 


m + r — fr 
m + r + fr 


k “in 


m + r 


(25) 


if  edge  effects  are  neglected.  This  is  a one-dimensional  linear  Bayes  filter  of  the  type  described 
in  the  preceding  section.  When  applied  successively  to  the  rows  and  then  the  columns  (or  vice 
versa)  of  a two-dimensional  pixel  array,  it  results  in  a two-dimensional  linear  filter  of  the 
exponential  window  type  described  in  Section  ll-A,  with 

l-kl/l 

h(i.  J)  “ constant  x <>  “ 


Hence,  this  type  of  two-dimensional  filter  can  be  implemented  recursively  for  an  N x N image 
by  2N  successive  applications  of  the  one-dimensional  filter  of  Eqs.  (19)  and  (22).  Such  an 
implementation  has  the  property  of  at  least  being  a Bayes  filter  for  the  corresponding  one- 
dimensional image  described  by  (14)  and  (18),  and  is  extremely  simple  computationally. 


Unfortunately,  as  noted  earlier,  this  particular  two-dimensional  filter  is  not  isotropic. 
Ideally,  one  would  like  to  construct  a one-dimensional  Bayes  filter  like  (19)  and  (22)  that  would 
give  a Gaussian  one-dimensional  weighting  function  away  from  the  end  points.  Such  a filter 
would  result  in  an  isotropic  "Gaussian  window"  in  two  dimensions  when  applied  successively  to 
the  rows  and  columns.  This  goal  does  not  seem  quite  attainable,  but  arbitrarily  close  approxi- 
mations to  it  can  be  achieved  at  the  expense  of  resorting  to  ever  higher-order  Markov  processes 
in  place  of  that  of  Eq.  (IS).  For  example,  the  next-higher  (second)  order  approximation  can  be 
constructed  by  defining  a Markov  process  with  state  vector  x such  that 
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and  by  giving  x the  dynamics  of  a "critically  damped  second-order  system"  driven  by  white 
noise.  The  recursions  analogous  to  (19)  and  (22)  in  the  resulting  Bayes  filter  for  the  x 
sequence  now  involve  2x2  matrix  operations.  The  form  of  the  equivalent  one-dimensional 
convolution  function  for  the  m’s  is  shown  in  Fig.  3 for  the  limiting  case  of  very  small  pixel 
spacing.  For  comparison,  this  figure  also  shows  the  Gaussian  window  and  exponential  window 
(the  first-order  approximation),  that  give  the  same  degree  of  grain-noise  suppression.  It  is 
apparent  that  this  second-order  approximation  is  indeed  closer  in  form  to  the  Gaussian  window 


Eig  3 — Convoluiion  functions  for  one-dimensional  fillers  with  equal 
grain-noise  suppression 


III.  A CLASS  OF  NONLINEAR  BAYESIAN  FILTERS 
A.  Formulation 

All  Bayesian  filters  described  in  detail  so  far  have  been  for  the  case  of  density- 
independent  grain  noise.  Although  the  additive,  Gaussian,  and  spatially  uncorrelated  character 
of  photographic  grain  noise  is  fairly  well  confirmed  by  experiment,  its  rms  magnitude  is  known 
[1]  to  increase  with  the  optical  density  of  the  original  image  as  approximately 
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where 


(r(D) 


o-(Do) 


D » film  density  at  the  point  in  question 
Do  — reference  density  (usually  1) 


a — rms  grain  noise  value. 


(26) 


In  this  section  we  seek  to  construct  Bayesian  filters  that  take  this  density  dependence  into 
account.  They  will  be  developed  according  to  a perturbation  analysis  for  small  D/Dq,  however, 
and  will  be  only  first-order  approximations  to  exact  Bayesian  filters.  The  exact  form  of  the  a 
priori  probability  distribution  that  should  be  used  for  the  original  image  is  a more  imponderable 
issue,  and  there  is  less  reason  to  depart  from  the  convenient  Gaussian  distribution  described 
earlier.  However,  there  is  some  motivation  for  dealing  with  spatial  image  fluctuations  that  are 
Gaussian  in  the  intensity  variable,  rather  than  density,  because  many  optical  processes  result  in 
image  transformations  that  are  linear  in  terms  of  intensity,  and  the  class  of  Gaussian  probability 
distributions  is  closed  under  linear  transformations.  Dealing  with  such  cases  leads  to  the  inclu- 
sion of  either  measurement  nonlinearities  or  nonlinear  recursion  relationships,  depending  on 
whether  the  analysis  is  performed  for  the  optical  density  or  intensity  variable. 

These  considerations  lead  to  the  formulation  of  a general  problem  in  which  a state  vector 
X evolves  sequentially  according  to  the  dynamics 

x,  + i — F,x,  + w,  + r'‘x,x,  + / — 0,  • • • , N — 1 (27) 

and  generates  the  noisy  observations 

2,  - H,x,  + dijxiX,  + y,',  I - 0,  • • • , N (28) 

where  the  following  conditions  hold. 

• The  components  of  the  three-way  matrices  T,,  n„  and  <I>,  are  approximately 
infinitesimal,  say  of  order  h.  All  other  quantities  are  of  order  unity,  including  F,*'- 

• Given  xq.  ' ' ' . x^y,  the  w,  and  v,  are  statistically  independent  zero-mean  Gaussian 
random  variables  such  that 

cov(w,)  - Q,  + I'ff  iX, 
and 

cov(v,)  “ R,  + 2Y'x, 

to  first  order  in  A.  'R,  and  Y,  are  of  order  h,  and  Q„  R„  and  R,"'  are  of  order  unity. 

• Qi  and  R,  are  symmetric  and  positive-semidefinite,  R,  is  positive-definite.  F,  • P/, 

n,  - n/.  'F,  - and  Y,  - Y/.  N <<  -j-  and  ||Fj|  < 1. 

n 
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The  objective  here  is  to  determine,  at  least  to  first  order  in  /i,  the  conditional  probability  distri- 
bution of  X,  (for  each  / individually)  given  the  observations  zq,  ■ ■ ■ , Zn  when  Xg  has  an  a priori 
probability  density  of  the  general  form  given  by  the  following  expression  to  first  order  in  h : 


- 


l+j(x  - x)^(V-'AV-')'V-'(x  - x)(x  - x) 


-y  (t-J)'^V-'(i-i) 


(2ir)"'^|V|'/^ 


(29) 


In  the  above,  V and  A are  symmetric,  n is  the  dimension  of  x,  A is  of  order  h,  and  x,  V,  and 
V"'  are  of  order  unity.  This  problem  is  treated  in  multivariate  form  to  accommodate  the  kinds 
of  recursive  realizations  of  image  fluctuations  described  in  the  preceding  section. 


The  problem  will  be  treated  sequentially  by  first  determining  p(x„/Z,)  for  successively 
greater  values  of  /,  where  Z,  - (zg,  • ■ • , zj,  and  then,  starting  with  p(xjy/Z/v),  finding 
p(xJZs)  for  successively  smaller  values  of  /.  It  turns  out  that,  to  first  order  in  h,  all  of  these 
probability  densities  are  of  the  form  of  (29).  In  general,  Eq.  (29)  can  assume  negative  values 
for  sufficiently  large  magnitudes  of  (x  — x)  and  must  be  modified  slightly  to  be  a proper  proba- 
bility density.  Because  of  the  rapid  decay  of  the  exponential  factor  and  the  orders  of  magnitude 
assumed  for  the  parameters,  however,  such  modifications  can  be  confined  to  a region  whose 
probability  mass  is  negligible  to  arbitrary  order  for  sufficiently  small  h.  Hence,  for  brevity,  Eq. 
(29)  itself  will  be  treated  as  a proper  density.  Since  it  has  the  form  of  a Gaussian  density  mul- 
tiplied by  a polynomial,  it  follows  from  standard  results  for  Gaussian  moments  that  its  integral 
over  R''is  unity  and  its  first  three  central  moments  are  given  by 

E(x)  - X -I- Tr(V-'A),  (30) 

cov(x)  — V,  (31 ) 


and 


E[(x  - E(x))(x  - E(x))^(x  - E(x))')  - 2 A 


(32) 


to  first  order  in  h.  Furthermore,  if  y — Ax  -I-  6,  where  A and  b are  constants  and  A"'  exists 
and  is  of  order  unity,  then  p(y)  is  also  a probability  density  of  form  (29),  with  the  parameters 
X,  V,  and  A replaced  respectively  by  Ax  -I-  b,  AVA^and  (AAAO'A^  It  is  also  useful  to  note 
that  the  apparently  more  general  probability  "density" 


p(x)  - 


1 -f  X ^(x  - a)  -f  y tr  {L((x  - a)(x  - a)  VJ)  -I-  y (x  - a)  '^Afx  - a)(x  - a) 


(27r)'’^^|V|''2 


(33) 


in  which  L is  symmetric  and  K and  L are  of  order  h,  can  be  reduced  to  the  form  of  (29)  to  first 
order  in  h.  This  is  done  by  making  the  substitutions 

X - a + VX,  (34) 


V - V -(-  VLV, 


(35) 


and 


A - (VAV)'V, 


(36) 


and  expanding  the  exponential  and  determinant  functions  in  Eq.  (33)  about  the  values  xand  V. 
Incidentally,  if  either  form  is  expanded  about  its  mean  and  covariance  to  first  order,  the  result 
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is  a standard  first-order  Edgeworth  expansion  of  p{.x).  Thus,  an  earlier  proposal  of  Sorenson 
and  Stubberud  |9l  to  approximate  posterior  density  functions  by  Edgeworth  expansions  in  Baye- 
sian filtering  is  accurate  to  first  order  in  this  context. 

Now,  if  > — X r"xx,  with  F — P^of  order  /i,  then  either  x — y — r"yy  to  first  order  or 
|x|  » 1 for  any  y of  order  unity.  The  Jacobian  matrix  (dx,/dy,l  can  be  expressed  as  I -f  2r"x, 
whose  determinant  is  1 + 2tr(rx)  to  first  order.  If  x has  the  density  (29),  then  the  standard 
formula  for  the  transformation  of  probability  densities  gives 

P(y)  - |l  + y (y  - x)  ^(V-'AV-')'V*'(y  - x)(y  - x)  - 2 tr  (Fy) 


to  first  order  for  y of  order  unity.  This  follows  because  large  values  of  x have  negligible  proba- 
bility, and  y — X is  otherwise  of  order  h.  Equation  (37)  is  also  an  adequate  approximation  for 
larger  y because  x cannot  then  be  of  order  unity.  Expanding  the  exponential  factor  of  Eq.  (37) 
about  (y  - x)  and  using  Eqs.  (33)  — (36),  we  note  that  p(y)  is  again  of  form  Eq.  (29)  to  first 
order,  with  the  parameters  x,  V,  and  A from  p(x)  changed  respectively  to 


X r 'xx  - 2VTr(F'), 

(38) 

V + 2(VF'  + F "V)'x, 

(39) 

A + VFV  + (VFV)'  (VFV)". 

(40) 

{ We  proceed  from  this  point  by  considering  two  arbitrary  adjacent  indices  / and  i + \ , 

dropping  for  brevity  the  subscripts  from  x„  F,  , Q„  w,,  and  fl,,  and  denoting  x,  + | by  y.  First, 
p(y/Z,4|)  is  derived  to  first  order  as  a probability  density  of  form  (29)  for  pix/Z,)  also  of  this 
form,  giving  a forward  recursion  on  /.  Then,  for  p(y/Z/v)  of  form  (29),  the  density  pixlZ^)  is 
derived  as  another  density  of  this  form,  which  specifies  a reverse  recursion.  Since  p(xo)  is  of 
form  (29),  these  two  results  together  establish  that  all  p(x,/Zv)  have  this  form,  and  they  pro- 
vide formulas  for  computing  the  parameters  of  these  densities.  It  is  possible  that  the  random 
variables  z,  can  by  chance  assume  such  excessively  large  values  that  the  assumed  approxima- 
tions are  not  valid,  but  this  only  happens  with  negligible  probability  for  sufficiently  small  h. 

B.  Forward  Recursion 

Suppose  that  the  probability  density  p{x/Z,)  is  of  form  (29)  with  parameters  x„  V„  and 
A„  and  with  corresponding  mean  denoted  by  x,.  Letting  t denote  Fx  implies  that  / has  a proba- 
' bility  density  of  the  form  of  Eq.  (29)  with 

; X - Fx„ 

I V-FV,F^ 

I 

A - (FA,F'')'Fr 

Letting  s - Fx  -I-  F''xx  with  F - F'^,F''  implies  that  j - / + F'  n.  Hence, _the  density  of  s is 
also  of  the  form  (29)  to  first  order,  with  corresponding  parameters  s,  V and  A given  by  (38)  — 
(40)  as 


(2ir)"''lV|'^^ 
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t 


J - X + ria  - 2VTr(r). 

V - V + 2(vr'  + r "v)'x, 

and 

A - A + vrv  + (vrv)'  + (vrv)" 

If  Q*'  exists,  it  follows  from  a first-order  expansion  of  the  exponential  and  determinant  in 
p(w/x)  that 


p(w/s)  - p„,^(w,F~'s)  - ll  - tr|0'''|r's(I  - Q^'h/wOH 


> * 


to  first  order  in  h,  where  'V  - For  /■  - w -t-  (r'ww.  Eqs.  (37) -(40)  imply  that 


p(r/s)  - 


1 - Tr^(n'  n")r  - tr(Q-''jr's(I  - Q-'r/-0) 


-I-  y rMO-'ft"  + n'Q-'  + (Q-'n")'lrr 

to  first  order.  Since  >-  — a -»■  j by  definition, 

p{y/Z,)  - - |■)Pr/^'^y  - r)dr 


3 * 


(2rr)"'2|Q|i/2 


k(r) 


- i|r''Cf'r  + ( y-r-K)^p-Uy-r~i)\ 

^ dr 


(41) 


(42) 


(27r)''|VQ|''^ 

to  first  order  in  h,  where  the  polynomial  factor  k(r)  is  given  by 

<t(A)  - 1 -K  (^-r-J)^(V-'AV*')'V-'(>'-r-s)(>'-r-s)  - Tr^(n'-I- ft")A 


(43) 


- tr(Q'''F'(>-A)(l-Q-'rrOl  + "-t-  a'Q-'-l- (Q-'fl  ")'Irr. 


Completing  the  square  in  the  exponent  of  the  integrand  gives 

(2ir)"^2|M|'^^  (2ir)"'^|V  - VM-'V|''^ 

where  M — V -(■  Q.  This  result  does  not  involve  Q*',  and  in  fact  holds  for  singular  Q as  well. 
In  the  latter  case,  however,  a more  complicated  derivation  is  required,  in  which  the  integration 
of  Eq.  (43)  and  the  expansion  corresponding  to  (42)  are  performed  over  the  subspace  of  R'for 
which  p(r/s)  is  nonsingular.  The  integral  of  (44)  has  the  form  of  the  expected  value  of  a 
third-degree  polynomial  in  r and  (>— D with  respect  to  a Gaussian  distribution  for  r whose 
mean  is  proportional  to  (y  — J).  Therefore,  it  can  be  expressed  in  the  form 

constant -I- \''(>'-J) -I- tr  (LKj; -?)(>' -J) Mil -h (j;  - s) '^(M A •M)’M(>- - Dfy  - s), 

where  the  constant  is  independent  of  >>  — s!  It  can  be  shown  that  L and  A ‘Mn  be  taken  as  sym- 
metric here  because  M is,  and  that  X,  L,  and  A * are  of  order  h because  A and  'F  are.  Since 


> 


I 


t 

( 


15 


I 


i 

i W W WILLMAN 


« 


p(ylZ)  is  a probability  density,  the  constant  term  must  be  unity.  Thus,  p(ylZ)  is  of  the  form 
(33),  and  therefore  reducible  to  the  form  of  (29)  to  first  order.  The  first  three  central 
moments,  which  determine  this  density,  can  be  found  by  performing  the  integration  of  Eq.  (44) 
and  substituting  the  result  into  the  earlier  definitions.  However,  it  is  simpler  just  to  evaluate 
them  by  using  (27)  and  decomposing  the  expectations  into  marginal  expectations  of  conditional 
expectations  given  x.  This  gives  the  result 

E(>-/Z,)  - Fx,  + Tr(r,(jf,V+  V,)  + HQ)  - x,.^  (45) 

co\(ylZ)  - FV,F^+  Q + 2('F,  + FV,r,  + r;v,FO'x,  (46) 

y E((>--x,^,)(>-x,^.,)^(>-x,.„)7Z,l  - (FA,F0'F^+  0 + 0'  + 0”; 

where 

0 - FV,r,V,F^+ 9';V,F''+ QHQ.  (47) 


Turning  now  to  the  updating  of  this  conditional  density  with  the  next  measurement 
we  delete  for  brevity  the  subscripts  of  v,+i,  Hj+j,  R,+i,  and  Y,+|.  Since,  to  first  order 
p(vly)  is  mean-zero  Gaussian  with  covariance  matrix  R + 2Y'y,  p(z/ y)  is  Gaussian  with  the 
same  covariance  matrix  and  mean  Hy  -t-  ^"yy.  Expanding  the  exponential  and  determinant  of 
this  density  function  to  first  order  about  the  values  indicated  below,  we  obtain 


p(z/y)  - (1  -I-  (r-Hy)^R-'«J)>>-  - tr  |R-'Y>[I -R-'fr  - Hy)(z -Hy)  Hll 

X i— 

(2rr)'"/2|R|''2 


where  m is  the  dimension  of  z.  Since  v is  independent  of  the  preceding  noise  values,  the 
updated  conditional  state  density  is  proportional  to  the  product  p(zl y)p(ylZ)  as  a function  of 
y,  which  follows  from  the  Bayes  rule.  Recalling  that  p(ylZ)  is  the  probability  density  specified 
by  (45) -(47),  completing  the  square  in  the  exponent  of  this  product,  and  using  the  "matrix 
inversion  lemma*  (6),  we  obtain 


p(y/Z,+|)  - g 


1 + y (y  -y)(V-'AV-')'V-'(y  -y)(y  -y) 


- triR-'Y'yll  - R'Kz  - Hy)(z  - Hy)^) 


e 


(2ir)"'^|V|''* 


where  g is  a constant  of  proportionality. 


(48) 


V - V - VH^(R  -I-  HVH0"'HV, 
a - y VH^R-'(z  - Hy), 

and  y,  V,  and  A are  now  used  to  denote  the  parameters  of  piy/Z,)  as  defined  by  (29).  The 
polynomial  factor  of  (48)  can  be  expressed  to  first  order  as 

l-l-X^(y-a)  + y tr|L((y  -a)(y  - a)*"- VI)  + y(y-a)^A(y-a)(y-a),  (49) 


J 


I 


i 

i 


I 
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in  which 

X - (2«|)  - H''R-'(Y  + 2HV4»')'1J;(R  + HVHO-‘(r  - Hy)  - H''R-Vjiv 
+ Tr{[H(V -•  A)'H 2HV<|) " + YJ  (R  + HVH  0 "'(z  - Hy)  (z  - Hy)  ^ 

X (R  + HVH  0 ■')  - H ^R-'Tr[HV<|>VH  ^(R  + HVH  0 ''(z  - H>)  (z  - HJ) 
x(R  + HVHO-'|-Tr(YR-'). 

i 

L - 2[(V-'AV-')'H <!»'  - H '■R-'(Y  + 4) "VH 0 + (H '"R-'YR-'H  - H '■R-'<t")'VH ^ 

- ( YR-‘H)  "1  (R  + HVH  0 -'(z  - HJ)  ‘ 

+ 2(H  ^R-'YR-'H  - H ^R-'<I>"  - ♦■R-'H)7, 

and 

A - (V-'AV-')'V-'  + (H^R-'Y  - <I>')R-'H  + [(H^R-'Y  - <D')R-'H1' 

+ [(H^R-'Y  - 4>')R-'H1". 

Consequently,  according  to  (33)-(36),  ;>(y/Z,+,)  can  be  reduced  to  the  form  of  (29).  If  we 
follow  the  suggestion  of  Sorenson  and  Stubberud  [101,  we  can  express  the  first  three  central 
moments  of  this  density,  denoted  here  as  V,+,  and  2A,+i,  to  first  order  as 

- x,+i  + V(H  + 2<I>"x,+,)  ^[R  + 2Y'3f,+,  + (H  + 2<I)'’x,+,)  V(H  + 2«l>"jf,+,)  H 
X (r  - Hx,+,  - 4)  "x,+|X,^.,  - Tr(V4>)]  + VTr|(R  + HVH  0 ‘'(Y  + 2HV(<t> 

- H '■R-'Y')  + H(V-' A)'H  n (R  + HVH  0 -'(z  - Hx,+,)  {z  - Hx,+,)  ^ 

- (R  + HVH  Oil  - VH  ^R-'  X Tr{(R  + HVH  0 ‘'HVitVH  ^(R  + HVH  O'' 

X ((z-Hx,+,)(z-Hx,+,)'^-(R  + HVH01I,  (50) 

V,^,  - V - V(H  + 2«l>''x,^,)  Or  + 2 Y'x,^, 

+ (H  + 2<I>"x,+,)  V(H  + 2«I)'’x,^,)  0 -'(H  + 2<I>’’x,+,)  V 

'I  + 2{(VV-'AV-'V)'H'^+(V(H^R-'YR-'H-4»'R-'H-H^R-'«|)")V1'VH^ 

+ (V(<J»-H'’R-'Y'-Y"R-'H)Vl')(R  + HVH0-'(z-H3f,+,),  (51) 

and 

A,+,  - (VAV)'V.  (52) 

This  result,  together  with  those  of  Eqs.  (45)~(47),  is  the  same  to  first  order  as  that  proposed 
I on  more  ad  hoc  grounds  by  Sorenson  and  Stubberud  [10]  for  obtaining  E(x,+,/Z,+,)  and 

I cov(x,+|/Z/+i)  from  E(x,/Z,)  and  cov(x,/Z,),  except  for  the  inclusion  of  the  measurement- 

driven  term  in  (51)  and  the  last  two  terms  of  (50)  containing  the  scatter  matrix  of  (z  — Hx,+i). 

The  use  of  such  scatter  matrix  terms  was  proposed  by  Lee  and  Tung  (11)  in  a closely  related 
I context. 

1 

il 
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C.  Reverse  Recursion 

We  now  seek  to  determine  p(xJZn)  as  a probability  density  of  the  form  of  (29)  to  first 
order,  given  that  p(xJZ)  and  p(x,^.,/Z/v)  are  both  of  this  form.  For  brevity,  the  latter  two 
densities  are  denoted  as 


and 


p(xlZ)  - 


pfy/Z/v)  - 


1 + y (x -x)  ^A(x  - jr)(x  - x) 


1 + -y)(y  -y) 


- i («-jr)'^V-'(*-«) 
> * 


(27r)'’'2|B|'^2 


(53) 


(54) 


The  results  for  the  general  case,  although  obtainable  by  the  methods  used  below,  are  very 
cumbersome,  and  their  interpretation  is  still  not  well  understood.  Hence,  the  development  is 
limited  from  this  point  to  the  case  of  zero  F,,  ft,,  and  Nk,,  which  simplifies  the  dynamics  of 
(27)  to 


y - Fx  + w;  w is  Gaussian  (0,Q) 


(55) 


with  the  omission,  as  before,  of  subscripts.  From  the  Markov  properly  of  the  x’s,  p{x/ y.Z^) 
- p(x/y,Z,),  so  that 


p(xlZs)  - / p(x/y.Z,)  p(ylZn)dy. 


(56) 


where  the  integration  is  performed  over  the  subspace  of  R"for  which  pfx/y.Z,)  is  nonsingular. 
For  brevity,  we  again  treat  only  the  case  of  nonsingular  Q,  in  which  case  this  subspace  is  R" 
itself  However,  the  result  does  not  involve  Q'',  and  can  be  verified  for  singular  Q as  well  by 
resorting  to  a more  complicated  derivation.  From  the  Bayes  rule,  p(x/y,Z,)  is  proportional  as 
a function  of  x to  the  product  pfy/x.Z,)  pix/Z).  Furthermore,  the  Markov  property  of  the 
X 's  can  be  exploited  again  to  replace  the  first  factor  in  this  product  by  p(y/x).  From  Eq.  (55), 

1 


piyix)  - 


( y-rxl^Q-‘<  y-ft) 

(27r)'/2''(Q('/i  ■ 


(57) 


The  product  of  functions  (53)  and  (57)  can  be  normalized  so  that  its  integral  is  unity  by 
completing  the  square  in  the  exponent,  which  leads  to 

p(x/y,Z,)  - jl  + (A  + Am  + Z"uu){x  - St  - VF^M^'m) 


+ tr 


T L + M'm 
2 


((x  - X - VF^M-'m)(x  - X - VF^M-'m)''-  si 


+ y (x-x-VF'^M-'M)'A(x-jf-VF'M-'M)(x-x-VF^M  '«) 


J((-jf-VF''M  /-VF'^M  'u> 

(2»r)''''^|S|''' 


(58) 
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where 

u ^ y -y. 

X - X + VF^M-‘(J  - Fx), 

M - FVF^+Q, 

S - V-VF^M-'FV, 

IT  - (M-'FVA)', 

I - M-'FVAVF^M-', 

X - Tr[Z(J-Fx)(7-Fx)n, 

L - 2AVF^M-‘(J  - Fx),  and 

A - 2r'(jr-Fx). 

Using  (54)  and  (58)  in  (56),  completing  the  square  in  the  resulting  exponential  factor, 

and  deleting  second-order  terms  leads  eventually  to  the  equation 


P.-aitlZy)  - 


e ^ 

f 

(2ir)''/2|G|'/2 

1 + K^t  + ytrlL(//^-  S)) 


+ -j-  t^Sti  - IM-’FVX  + Tr(Sn)]^u  - y tr(M’'FVAMMO 

-I-  y M^(e  - rvF^M -')uJ 

where 

G - V - VF^M-'(M  - B)M-'FV 
and 

K - B - BM-'FVG"'VF^M''B 

- Q(MB"'Q  -I-  FVF0“'M  (by  the  matrix-inversion  lemma). 


e ^ 

(2ir)"/^|K|'^^ 


(59) 


(60) 

(61) 


The  integrand  of  (59),  as  a function  of  u,  has  the  form  of  a third  degree  polynomial  multiplied 
by  a Gaussian  density.  Therefore,  the  integral  can  be  evaluated  by  using  standard  results  for 
Gaussian  moments.  This  procedure  gives 
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p(xlZs)  - 


1 + \^(x  - x)  + y tr(L|{x  - x)(x  — x)  ^ — Gl) 


+ y (x  - jf)  ^A(x  - x)(x  - x) 


(27r9"^^|G|'^^ 


where 


X - X + F^O-'K[Tr(Ke  - KE'VF^M-')  - M-'FV\  - Tr(Sn)l. 


(62) 

(63) 


L - L - F'Q-'KM-'FVLVF^M-'KQ-'F  (64) 

A - A + F'Q-'K(0  - rVF^M-')KQ-'Fl'KQ-'F.  (65) 

Since  p(x/Z  v)  is  of  the  form  (33),  it  is  reducible  to  the  form  of  (29)  to  first  order.  It  is  con- 
venient at  this  point  to  define  the  more  meaningful  variables. 


A — (VAV)'V  (1/2  third  central  moment  of  pix/Z,)), 

H — (B(9B)'B  (1/2  third  central  moment  of  pf^/Z/v)), 

X - x-(-Tr(V''A)  (mean  of  p(x/Z,)), 

y “ >-l-Tr(B''S)  (mean  of  pfj'/Zv)). 


The  application  of  (33)  — (36)  to  (62)  then  shows,  after  much  manipulation,  that  the  first  three 
central  moments  of  p(x/Z/v)  can  be  expressed  to  first  order  as 


X + VF''M"'(.V  - Fx)  + G[V-'  - F^(MB-'Q  + FVF0’‘F] 

X Tr{M"'FAF’"M"'[(;!  - Fx)(>  - Fx)  (M  -t-  B)]),  (66) 

G -I-  2(GV-'AV-'G  - GF^(MB-'Q  FVF0-'FAF^(QB-'M  -1-  FVFO-'FGI’ 

X F^M''(>>  - Fx),  (67) 

and 


2((kv-'av-'k)'v-'k-i-1kf^mb-'q-i-fvf0-'((mb-'sb-'m)'b-'m 

- (FAFO'Fn(OB-'M  -I-  FVF0'‘FK)'FK1.  (68) 


D.  Overall  Result 

It  is  convenient  for  further  perturbation  analyses  to  express  the  preceding  results  in  terms 
of  the  "nominal"  covariance  matrices  of  p(xJZ)  and  pixJZ/^)  for  the  case  in  which  the  pertur- 
bation parameters  are  all  zero.  If  these  matrices  are  denoted  respectively  by  P,  and  B„  it  fol- 
lows from  standard  Kalman  filtering  and  Bayesian  smoothing  results  [6]  that  they  are  generated 
to  first  order  by  the  recursions 
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M,+i  - F,P,F,’^+ Q,;  Mo  specified  a priori,  (69) 

P,  - M,  - M,H,^(H,M,H,^+  R,)-'H,M,-2M,Tr(Rr'Y,Rr'Y,)M„  (70) 

and 

B,  - P,  - P,F/Mr;,(M,^,  - B,^,)Mr;,F,P,;  B^  - P,v-  (71) 

Let  us  define  the  following  quantities: 

D,  - y (cov(x,/Z,)  - P,1  (72) 

E,  - y lcov(x,/Zv)  - B,1  (73) 

A,  - y third  central  moment  of  p(xJZ)  (74) 

0,  — y third  central  moment  of  p(x,/Zv)  (75) 

K,^,  - B,^,  - B,^,Mr.',F,P,B,-'P,F,'’Mr;,B,^,  (76) 

X,  - E(x,/Z,)  (77) 

X,  - E(x,/Z^)  (78) 


and  suppress  for  brevity  in  the  following  notation  the  subscript  / for  the  variables  x,  F,  F,  P,  Q, 
0,4',  n,  and  A,  and  the  subscript  /+1  for  the  variables  M,  H,  R,  <I>,  Y,  0,  E,  B,  K,  x,  and  z. 
Then,  for  the  general  case  with  nonzero  perturbation  parameters,  the  preceding  recursions  can 
be  expressed  to  first  order  as  follows: 

Forward  Recursion  (/  - -1,  0,  • • • , N) 

X - Fx  + Tr[r(xx^+  P)  + flQl;  / > 0 
D - FDF^+ ('P  + FPF' + F'PFO'x;  / >0 
A-(FAFO'F^+FPrPF^+'P'PF^+QnQ  + (FPrPF^+'P'PF'^+OnQ)' 

+ (FPrPF^+  'P'PF^+  QftQ) / > 0 

x,+,  - X + P,+|H  ^R-'lz  - Hx  - <I>''jcx  -Tr(M<l>))  + 2P,+,(M-'DH ''+  <l>x 

-H^R-'(HM<I>'  + Y)'x](R  + HMHO''(z-Hx)+P,+,Tr((R  + HMHO'' 

X III  + (YR + R-'Y)'x]  Y + 2HM(<I>  - H '^R-'Y')  - (H  ^R-'(HM<t>MH  H")' 

+ H(M-'A)'H  n(R  + HMH  0''((r  - Hx)(z  - Hx)  (R  + 2Y'x  + HMH  011  (82) 

D;+,  - P,^,M-'DM-'P,+,  + IP,+,(H  ^R-'YR-'H  - <1>'R-'H  - H '^R-'«l>  ")P/+il'x 
+ |(P,+,M-'AM-'P,+,)'H'' 

+ [P,+,(H^R-'YR-'H-<I>'R-'H-H^R-'<D")P,+irMH^ 

+ IP, +,(«!>  - H ^R-'Y'  + Y "R-'H)P,+iri(R  + HMH  0 "'(r  - Hx)  (83) 


(79) 

(80) 

(81) 
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A,^,  - (P,^,M-'AM-'P,^,)'M-'P,+,  + (P,^,H''R-'YR-'HP,^,)'P,^, 

-(P,^,<t.P,^,)'R-'HP,^, 

+ 1(P,^,H''R-'YR-'HP,^,)'P,^,-(P,^,<I>P,^,)'R-'HP,^,1' 

+ [(P,^,H^R-'YR-‘HP,+,)'P,+,  - (P,^,<l>P,^,)’R-'HP,+,l"  . (84) 

Reverse  Recursion  for  F,  'P,  (1  — 0 (/  — N — 1,  •••,0) 

X,  - x + PF^M-'(x-Fjf)+2(I-PF^M-'F)DF'M-'(x-Fjf) 

+ B,[p-'-F^(MB-'Q  + FPFO-'F]TrlM-'FAF^M-'((x-Fx)(x-Fx)'' 

-(M+B)]l;  jfA/-Xv  (85) 

E,  - D-(1-PF^M-'F)DF^M-'(M-B)M-'FP-PF^M-'(M-B)M-'FD(I-F'^M-'FP) 

+ PF^M-'(FDF^-  E)M-'FP  + (B,P-'  AP-'B,-B,F^(MB-'Q  + FPF')'' 

xFAF^(QB*‘M  + FPFO'‘FB,]F’'M-'(x-Fx);  E^  - (86) 

e,  - (KP-'AP-'K)'P-'K  + |KF'(MB-'Q  + FPF^)-‘l(MB-'eB-'M)'B-'M  > 

-(FAF'')'Fn(QB-'M-FPFO''FKrFK;  Wv-Av.  (87) 

For  the  index  / — -1,  the  quantities  x,  D and  A are  specified  by  the  a priori  distribution  of  xq, 
not  by  (79)  — (81).  For  the  common  limiting  case  of  a "fiat  prior"  (i.e.,  an  a priori  probability 
density  for  Xq  which  has  an  approximately  infinite  diagonal  covariance  matrix),  a reasonable 
procedure  for  initiating  the  recursions  here_would_be  to  make  Mo  a diagonal  matrix  with  very 
large  diagonal  entries,  and  use  zero  for  x,  D and  A in  the  initial  iteration  of  (82)  — (84).  For 
the  special  case  of  Hq  ” I,  this  procedure  could  be  circumvented  entirely  by  using 

Pfl  “ R 0-  ( 

Xo  “ Zfl  ~ *^*0  ^0^0  “ RoTr(YoRo  '). 

Do  ” (Yo  ~ Rn‘J*o  ~ ‘I’o  Ro)zo.  and 

Ao  — Y oRo~  Ro‘l’oRo+  (Y oRo“  Rn*^oRo)’  + (Y o'Ro”  Rn'l’oRo)  "■ 

Since  all  posterior  probability  densities  p(x,/Z^)  are  of  the  form  of  (29),  they  are  com- 
pletely specified  to  first  order  by  their  first  three  central  moments,  which  in  turn  are  specified  to 
first  order  as  x„  B,  + 2E„  and  2W,  by  the  system  of  recursions  (69)  — (71)  and  (79)  — (87). 

Asymptotic  Results 

In  practice,  the  pixel  spacing  is  often  so  small  that  the  number  of  pixels  is  no  longer  small 
compared  to  \/h  for  the  nonlinearities  that  are  present,  in  such  cases,  one  is  typically  faced 
with  a situation  in  which  the  problem  parameters  assume  the  form 

F,  - I -t-  «F(/) 
r,  - fl'lr) 

ii,  - n(/) 
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H,  - H(0 

<l>,  - <I>(0 

R,  - - R(/) 

€ 

Q,  - «Q(/) 

Y,  - - Y(/) 

< 

where  « denotes  the  small  quantity  1/N  and  r denotes  the  quantity  //N,  a variable  ranging  from 
0 to  1.  As  long  as  « > > a consistent  derivation  of  the  above  results  can  still  be  developed, 
in  this  case  by  folllowing  the  sort  of  procedure  used  in  Refs.  12  and  13,  except  that  all  terms  of 
order  /i€  or  larger  are  now  retained.  This  analysis  becomes  extremely  tedious,  however,  and  is 
not  exhibited  here. 

IV,  EXPERIMENTAL  RESULTS 
A.  Filtering  in  Density  Variable 

To  evaluate  the  performance  of  these  filters  experimentally,  we  chose  a test  image  as 
similar  as  possible  to  one  used  by  Naderi  [1]  in  work  with  some  related  filters.  Our  image  was 
generated  by  scanning  a photographic  negative  on  a comparatively  grainy  film  (Kodak  Tri-X) 
with  a microdensitometer,  using  a S-^m  x S-^m  aperture  to  produce  the  128-pixel  by  128-pixel 
digital  image  shown  in  Fig.  4.  The  optical  densities  of  the  light,  medium,  and  dark  areas  of  this 
image  average  2.1,  1.2,  and  0.3,  respectively.  The  appearance  and  statistical  parameters  of  this 
example  agree  fairly  closely  with  those  of  Naderi,  which  were  entirely  computer-generated,  in 
addition  to  the  photographic  grain  noise,  our  image  undoubtedly  contains  some  shot  noise  con- 
tributed by  the  microdensitometer,  but  this  is  believed  to  be  relatively  minor  and  is  not  con- 
sidered further. 

Further  following  Naderi’s  approach,  we  constructed  our  filter  by  combining  one- 
dimensional filters  for  the  pixel  rows  and  columns,  each  of  which  is  treated  as  a noisily 
observed  Markov  process  vector  with  first-order  correlation  of  0.85.  This  means  that  the  origi- 
nal one-dimensional  "images"  can  all  be  regarded  as  generated  by  a simple  scalar  process  of  the 
form  of  (27), 

- /jf,  w,;  w,  is  Gaussian  (0,^)  (88) 

provided  that  the  Markov  process  is  assumed  to  be  stationary.  In  the  above, 

X,  “ u,  — u, 

u,  — optical  density  of  the  /th  pixel, 

u — average  pixel  density  of  the  entire  (two-dimensional)  image, 

/ - 0.85, 

q — 0-^(1  — / ^)  and  1 

I from 

<t‘  - variance  of  optical  density  fluctuations  in  [ Eq.  (17) 
the  two-dimensional  original  image.  J 


23 


W W WILLMAN 


T 


1 


In  addition  to  grain  noise,  optical  and  chemical  blurring  are  present  in  the  actual  photo- 
graphically recorded  image.  Naderi  estimated  the  magnitude  of  this  blurring  and  took  it  into 
account,  but  the  resulting  contribution  to  the  final  filter  design  turned  out  to  be  relatively  minor 
in  this  case  because  the  correlation  coefficient  / between  adjacent  pixels  is  so  close  to  unity  that 
it  overwhelms  the  compensation  in  the  design  for  the  expected  optical  and  chemical  blurring. 
This  blurring  still  affects  the  filter’s  output  appreciably,  but  mostly  because  it  affects  the  input. 
Thus,  although  for  simplicity  these  blurring  phenomena  are  not  considered  in  the  present  filter 
design,  the  situation  in  this  respect  is  still  comparable  to  that  of  Naderi. 

The  grain-noise  model  adopted  here  is  such  that  the  variance  is  a linear  fit  to  that  of  (26) 
in  the  range  of  the  optical  density  levels  of  interest.  Hence,  the  observed  image  corresponding 
to  (88),  with  the  average  optical  density  u subtracted  for  convenience,  can  be  treated  as  a 
sequence  (z, ) such  that 

2,  — X,  + V,  (89) 

where,  given  the  x,  sequence,  the  v,  are  independent  zero-mean  Gaussian  random  variables 
with 

E(v,Vx,)  - r -I-  2t)X,  — r + 2t}(i/,  — i7)  (90) 

in  ihe  x,  range  of  interest.  For  the  image  examined  here,  r and  tj  were  actually  determined  by 
matching  the  values  given  by  (90)  and  those  observed  at  m,  — u ± <t,  for  which 

u - 1.2 


and 


tr  — 0.6. 


The  resulting  parameters  are  approximately 

r - 0.1 


and 

T)  0.04. 

Since  t)  is  so  small,  it  can  be  regarded  as  an  approximately  infinitesimal  parameter,  so  (88)  and 
(89)  represent  a simple  scalar  instance  of  the  class  of  observation  systems  analyzed  in  the 
preceding  section.  Hence,  the  Bayes  filter  for  this  system  can  be  found  to  first  order  as  a spe- 
cial case  of  the  recursions  (79)  — (87),  for  which  the  desired  output  is  the  x,  sequence  of  poste- 
rior mean  values.  Once  these  are  found,  the  average  optical  density  li  has  to  be  added  back  to 
the  result  for  each  pixel.  Neglecting  edge  effects  allows  the  M„  P„  B„  and  H,  variables  to 
be  treated  as  constants,  denoted  here  by  m,  p,  h,  X,  and  f>.  This  reduces  the  pertinent  equations 
to  the  following. 


- / V,  -1- 

£+2 

r 

m r 

A 

. M. 

r — ini 

4' 

+ 2 

ni  r + 2t)X, 

f'x, 


(z,n  - f.i,) 


Kz, H - / iy  - (w  + r + 2t}X,)1 


(91) 
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d(+\  “ 

d,  + V 

4/ii.  + 

2t) 

£ 

m 

"1 

r 

r 

(A+I  - fx)  (92) 


[(x,+|-/x,)^-(/M  +6)1.  (93) 

where  (91)  and  (92)  are  evaluated  successively  for  / - 0,  1,  • • • , N - 1 (where  N - 127  in 
the  present  case),  beginning  with 


X,  - X,  + 


X 

m 


P + 2 


Qd, 

m 


and 


•^0  - Zo 


do  - -nzo. 

and  (93)  is  evaluated  in  reverse  order  for  the  x„  beginning  with  X/y  — X/y.  The  constants  p,  m, 
b,  and  X are  determined  from  (69)  — (71)  and  (84)  by  requiring  that  the  variables  in  these 
equations  don’t  change  between  adjacent  indices.  This  yields 

P = [XR  + q-R  + + 9 - R)  + 4 R 

2r  J r + 2r)^ 

w ” f^P  + (f. 

^ _ p - i f p/m)^m 

1 - (/p/m)2  ’ 

and 

^ 3p\ 

rH\  - (/p/m)’l  ■ 

Note  that  the  parameter  0 is  not  needed  for  the  limited  purpose  of  computing  only  the  x, 
values. 

The  one-dimensional  Bayesian  filter  determined  by  recursions  (91)  — (93)  was  used  to 
construct  a two-dimensional  filter  for  the  image  at  hand  by  applying  it  successively  to  the  pixel 
rows  and  columns  in  the  manner  described  in  Section  IIC.  The  rationale  for  this  construction  is 
that  this  filter  is  a perturbation  of  the  linear  one-dimensional  filter  of  Eqs.  (19)-(23).  It  has 
already  been  shown  that  when  this  unperturbed  filter  is  used  on  a two-dimensional  image  in 
this  way,  the  result  can,  except  for  edge  effects,  be  interpreted  as  a two-dimensional  weighting 
function  of  the  reasonable  "exponential  window"  form.  Hence,  when  the  perturbations  are 
added  to  the  one-dimensional  filters  to  account  for  the  grain-noise  density-dependence,  it  is 
reasonable  to  combine  them  in  the  same  way.  It  should  be  noted  that  this  is  not  the  procedure 
used  by  Naderi  [1),  so  the  treatments  become  distinct  at  this  point. 

Figure  4 shows  the  result  of  applying  this  composite  filter  to  the  test  image,  which  is 
shown  in  the  same  figure.  For  comparison,  the  output  of  the  corresponding  unperturbed 
(linear)  filter  is  also  displayed.  The  constituent  one-dimensional  linear  filters  used  here  can  be 
obtained  simply  by  setting  t)  - rfo  - 0 in  (91)-(97);  therefore,  by  (90),  these  are  the  Bayes 
filters  resulting  from  approximating  the  grain  noise  as  independent  of  the  optical  density. 


(94) 

(95) 

(96) 

(97) 
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c.  Linearly  filtered  image 
Figure  4 
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Hence,  the  difference  between  Figs.  4b  and  4c  is  that  which  results  from  taking  the  density- 
dependence  of  the  grain  noise  into  account,  at  least  to  within  the  accuracy  of  the  linear  variance 
fit  used  here.  This  difference  is  not  dramatic,  but  it  is  noticeable.  As  might  be  expected  on 
intuitive  grounds,  the  more  accurate  representation  of  the  grain  noise  leads  to  relatively  more 
filtering  in  regions  of  high  optical  density,  where  the  noise  magnitude  is  larger.  The  range  of 
optical  densities  is  quite  high  in  this  example,  which  indicates  a good  degree  of  robustness  in 
the  approximations  and  perturbation  procedures  used  in  this  context. 

Direct  grain-noise  filtering  results  appear  unimpressive  partly  because  the  same  visual  per- 
ception tends  to  be  created  automatically  by  the  observer.  For  this  reason,  we  investigated  the 
value  of  such  filtering  as  a preliminary  operation  used  to  improve  the  subsequent  results  from 
more  intuitive  nonlinear  image  processing  procedures.  The  so-called  noise-cheating  algorithm 
described  in  Naderi  II]  was  used  for  this  purpose.  This  algorithm  is  essentially  a variant  of 
gray-level  slicing  with  an  extra  step  added  to  clean  up  boundaries  between  regions  of  distinct 
gray  levels.  Also,  the  gray  levels  to  which  the  pixels  of  the  input  image  are  quantized  are 
spaced  40-  apart,  where  o-  is  the  rms  noise  value  (possibly  a slowly  varying  function  of  gray 
level)  of  the  input  image. 

Figure  5 shows  the  results  of  applying  this  noise-cheating  algorithm  to  (i)  the  unfiltered 
test  image,  (ii)  the  linearly  filtered  image,  and  (iii)  the  nonlinearly  filtered  image.  In  these 
cases,  the  4o-  spacing  for  the  quantization  levels  were  determined  as  follows; 

(i)  By  direct  use  of  (26),  with  m — 1.2  and  a-(u)  yfr  0.31  (result  varies  with  gray 
level  u) 

(ii)  In  this  case,  the  result  of  (i)  was  multiplied  by  the  noise-suppression  factor  described 
in  Section  IIA  for  the  exponential  window  weighting  function  corresponding  to  this 
filter,  which  is  linear. 

(iii)  In  this  case,  the  action  of  the  nonlinear  filter  at  any  given  gray  level  u was  approxi- 
mated as  that  of  the  linear  filter  of  (ii)  when  the  grain-noise  variance  parameter  r is 
replaced  by  /•  -I-  2t)U.  The  computation  then  proceeded  as  in  (ii). 

The  appearance  of  the  end  result  is  affected  by  shifts  of  the  quantization  levels  in  each  case,  so 
several  different  choices  are  shown  for  each.  Although  no  general  improvement  in  appearance 
is  apparent  from  the  filtering,  it  should  be  remembered  that  the  quantized  gray  levels  in  the 
output  are  more  closely  spaced  in  the  case  of  the  filtered  images,  and  more  closely  spaced  in 
the  bright  areas  (and  less  so  for  the  darker  areas)  in  the  nonlinearly  filtered  image  than  in  the 
linearly  filtered  one.  For  example,  part  (iv)  of  Fig.  5 shows  the  effect  of  applying  this  algo- 
rithm to  the  unfiltered  test  image  with  the  more  closely  spaced  quantization  levels  that  were 
used  for  the  nonlinearly  filtered  image  in  part  (iii).  It  is  apparent  that  the  filtering  allows  a finer 
degree  of  gray  scale  resolution  for  roughly  the  same  degree  of  geometric  degradation  of  the 
image. 

B.  Filtering  In  Intensity  Variable 

For  reasons  cited  earlier,  the  estimated  image  is  often  desired  in  terms  of  the  intensity 
variable,  which  is  related  nonlinearly  to  the  optical  density  by  the  relation 
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tv.  UnHIlercd  test  image  — quanlizalion  tevels  of  part  ill 


Fig.  S — Applications  of  noise-cheating  algorithm 
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y - -logio  (p)  In  X - -0.434  In  x (98) 

where 

y • optical  density 

X - transmitted  intensity /incident  intensity  (—  1 here). 

Two  basic  approaches  to  this  problem  are  examined  here,  each  of  which  is  based  on  one- 
dimensional line  Alters  constructed  according  to  the  above  perturbation  analysis.  In  the  first 
approach,  the  observed  optical  density  is  regarded  as  a nonlinear  measurement  of  the  intensity, 
and  the  intensity  variable  is  filtered  directly.  In  the  second,  the  density  variable  is  filtered  after 
the  dynamics  of  the  intensity  fluctuations,  which  are  presumed  linear  and  Gaussian,  are  con- 
verted into  equivalent  nonlinear  dynamics  for  the  density  fluctuations.  In  both  cases,  the  opti- 
cal density  is  used  as  the  measurement  variable.  It  is  the  variable  that  behaves  experimentally 
as  one  corrupted  by  noise  that  is  approximately  additive  and  Gaussian,  as  required  by  the  per- 
turbation theory  used  in  the  filter  construction.  The  resulting  one-dimensional  filters  are  "cas- 
caded* as  before  to  produce  the  final  two-dimensional  filter.  The  work  here  is  fairly  preliminary 
and  is  limited  to  the  case  of  intensity  fluctuations  of  the  simple  form  described  by  Eq.  (88). 

1.  Nonlinear  Measurement  Approach 

In  this  approach,  the  dynamics  of  the  image  fluctuations  are  still  of  the  form  (88),  except 
that  the  u,  variable  now  denotes  (relative)  intensity  instead  of  optical  density,  and  u denotes 
average  intensity.  Subtracting  the  average  value  of  the  observed  optical  density  given  by  (98) 
as 

-0.434  jin  u, 

and  rescaling  the  difference  (and  the  associated  noise  variable)  leads  to  a measurement  variable 
of  the  form 

2,  - X,  + (/)x,^  + v„  (99) 

where  <p  is  chosen  to  approximate  the  logarithmic  nonlinearity  in  the  range  of  interest.  This  is 
again  an  instance  of  the  class  of  observation  systems  treated  in  the  last  section,  and  is  the  same 
as  the  preceding  case  of  optical  density  filtering  except  for  the  added  complication  of  a measure- 
ment nonlinearity.  The  filtering  equations  corresponding  to  (91)  — (97)  become  the  following 
in  the  present  ca.se: 

x,,,-/x, -i-  ^ U,,, - ./X,  - «(/V+  m)l 
(t)  + (w  - r)<^)x,j(z,  + |- /x,) 

fft  f > 

+ ^ + 2<fK/>/r)M(z,>i  - /x,)^  - (m  + r -f  2i7,r,)l.(l()0)  ; 

r r‘  m + r + 2t)X,  j 

tf,*!  “ ( fp/m)^d,  + (ri  - 2r<t>)(p/r)^fx,  i 

+ — - 2{ri- r<l>){p/r)^ ^ — (z,^i  - /x,).  (101)  ^ 

r w + r ? 
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and 

_ ipH-n  - r<l>) 

'■Ml  - (fp/m)^] 

The  equations  for  x„  p,  m,  and  b remain  unchanged  as  (93)  and  (94)  — (96). 

2.  Nonlinear  Dynamics  Approach 


(102) 


In  this  approach,  the  results  of  the  intensity  variable  hltering  are  expressed  in  terms  of  the 
optical  density  variable.  The  final  result  can,  of  course,  be  transformed  back  to  an  intensity 
representation  by  using  the  inverse  of  (98).  With  intensity  fluctuations  of  the  form  (88),  now 
expressed  as 

M,*^i  - u - /(«,  - m)  + (103) 

in  which 

u — average  intensity 

u,  - intensity  of  / th  pixel 

w,  — constant-variance  Gaussian  noise  sequence, 

a careful  expansion  to  second  order  in  the  image  fluctuations  shows  that  the  corresponding  opti- 
cal density  has  dynamics  given  by 

^,  + 1 - fx,  + «o,  - x}  + x,oi,  -I-  a>,M  (104) 

2a  a 2a 

where 

a - logiflf  ~ 0.434 

to  — - ^ w,  (another  Gaussian  noise  sequence,  variance  denoted  by  q) 
u 

X,  — / th  optical  density  — x,  x — — a In  u. 


Again,  the  zeroth-order  parameters  f and  q were  chosen  here  to  match  the  observed 
image  characteristics  (rms  density  fluctuations  and  degree  of  correlation  between  adjacent  pix- 
els). To  first  order  (104)  is  equivalent  to 

x,+,  - fx,  + ot,-  x^  -I-  m}  (105) 


with 

E(w,Vx,)  - q + 2 ^ X,. 

a 


(106) 


The  forward  recursion  of  (79)  — (84)  includes  the  present  case,  with  the  dynamics  of 
(105)  and  (106),  and  the  measurements  of  (89)  and  (90).  Neglecting  edge  efifects  by  using  the 
steady-state  values  of  p,  m.  and  q from  (94)  — (96),  and 


3p^ 

m’-/V  r, 

1 - (/p/m)’ 

2m’a  r\ 

(107) 
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leads  to  the  following  first-order  forward  recursion. 


■^,+1  - fx,  - (x^  f ~ 


-I-  2 


Jl. 


m r 

r 


3 

1 

1 

X, 

m 

ma  r 

(z,  + l - fx,) 


f\  + ^ (r-3w) 

2«  p(r  27)X,) 


(m+r)^ 

X [(2,  + |-/x,)^-(m -l-r  -I- 2t)X,)] 


^<4-1 


II 

2 

d,-¥ 

V 

£ 

2 

m -fp 

£_ 

fx,  + 

\2 

£ 

m 

r 

a 

m 

r 

r 1 

(2,  + l -/x,). 


(108) 

(109) 


The  general  reverse  recursion  described  by  (85)  — (87)  must  be  extended  to  generate  the 
x,’s  for  the  present  case.  There  is  no  conceptual  difficulty  in  simply  including  the  extra  terms 
in  the  analysis  there,  but  the  general  result  becomes  extremely  complicated.  Currently  we  are 
unable  to  simplify  it  due  to  insufficient  insight.  For  the  simple  special  case  at  hand,  however, 
these  details  have  been  worked  out  for  the  limited  purpose  of  computing  only  the  posterior 
means  x,.  The  result  is  the  following  reverse  recursion: 


jf X,  -I-  ^ 
m 


■ (/m  -(-  2fp  —3m)x, 


+ 2 

A. 

5 

d. 

fqK-\--^  if  pm  - m^-  3fpq ) 

m^ 

m^ 

2a 

X ((x,+i-/x,)^-(m  +b)]+  {/pm  - m^-3fpg)  + jf/v  “ x^-  (110) 

- ^ ma 


am" 


3.  Numerical  Reaults 

Both  of  these  approaches  were  applied  to  the  test  pattern  used  previously.  Unfortunately, 
the  logarithmic  nonlinearity  here  is  so  severe  that,  with  the  rather  extreme  contrasts  present  in 
this  particular  image,  the  first  procedure  would  not  work  without  some  serious  compromises 
being  made. 


In  the  first  procedure,  the  observed  optical  density  is  treated  as  a noisy  nonlinear  meas- 
urement of  the  intensity.  It  can  be  shown  from  (98)  and  (89)  after  much  manipulation  that 
(suppressing  the  / subscripts) 


z 


c -I- 


1 

0.868 


+ V, 


(111) 
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where 

X — pixel  intensity 

jf  - average  intensity  over  image  (—  0.063  here) 

7 - (y  - J)  + V 

> - pixel  optical  density 
y — —0.434  in  X, 


to  second  order  in  x - x,  the  departure  of  the  intensity  from  its  average  value  By  definition, 
the  intensity  value  can  be  recovered  from  the  f variable  as 


X 


X 


0.434 


(112) 


Ideally,  one  would  wish  to  use  - 1/0  868  in  the  perturbation  equations  ( 100) -(102)  for  filter- 
ing the  c variable  in  this  case,  so  that  (99)  and  (111)  agree  for  small  fluctuations  of  * How- 
ever, these  fluctuations  reach  values  of  ±0.8,  which  reduces  (112)  to  an  absurdity,  giving  pos- 
sibly negative  intensity  values.  Experimentally,  the  use  of  this  "low  contrast’  value  for  <6  led  to 
numerical  instability  when  used  on  the  test  image.  In  order  to  achieve  reasonable  experimental 
results,  we  found  it  necessary  to  reduce  <6  from  this  low-contrast  value  by  a factor  of  about  ten 
This  in  effect  sacrificed  low-contrast  accuracy  for  the  sake  of  a better  overall  fit  of  the  largo- 
scale  intensity  fluctuations  encountered  in  this  image  So  much  had  to  be  sacrificed,  however, 
that  the  use  of  this  tt>  parameter  to  represent  the  nonlinearities  here  is  of  doubtful  validity  and 
should  probably  be  confined  to  images  of  much  lower  contrast  The  output  of  the  filter, 
modified  as  just  described,  is  shown  in  Fig  6.  Rather  than  displaying  the  actual  intensity  vari- 
able specified  by  (112),  Fig.  6a  shows  the  filtered  values  of  c,  the  variable  on  which  the  filter 
actually  operates.  Since  (112)  is  linear  and  i is  the  conditional  mean  of  «,  the  conditional  mean 
of  the  intensity,  x,  is  given  by 


Hence,  this  result  is  the  same  as  the  negative  of  the  filtered  intensity,  at  least  to  within  an  arbi- 
trary bias  and  scale  factor.  Figure  6b  shows  the  result  of  applying  the  appropriate  linear 
transformation  to  the  i estimate  for  obtaining  the  filtered  image  in  the  origtnal  intensity  vari- 
able. The  exact  details  of  this  transformation  depend  on  the  dynamic  range  of  the  particular 
display  device,  etc.  Also,  an  optical  enlargement  of  the  original  photographic  negative  (i  e.,  the 
unfiltered  intensity)  is  shown  in  Fig.  6c  for  comparison. 


The  nonlinearities  here  did  not  have  such  a catastrophic  effect  on  the  second  method,  in 
which  the  filtering  is  performed  on  the  optical  density  variable  with  nonlinear  (one- 
dimensional) fluctuation  dynamics.  Here,  acceptable  performance  could  be  achieved  with  the 
nonlinearity  coefficients  at  roughly  their  low-contrast  values,  although  this  appeared  to  be  near 
the  borderline  of  feasibility  for  this  image.  Thus  it  appears  that  this  nonlinear  dynamics 
approach,  while  more  complicated,  is  basically  more  robust  than  the  nonlinear  measurements 
approach.  Figure  7a  shows  the  direct  result  of  using  this  filter  on  the  test  image,  which  appears 
in  terms  of  the  optical  density  variable.  Figure  7b  shows  the  same  result  after  having  been 
transformed  to  the  intensity  representation  by  applying  the  inverse  of  (98).  The  optical 
enlargement  of  the  original  photograph  is  repeated  for  convenience  in  Fig.  7c. 
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In  transforming  these  last  results  to  the  intensity  variable,  it  would  have  been  more  satis- 
fying to  compute  the  conditional  mean  of  the  intensity  itself  instead  of  the  nonlinear  (exponen- 
tial) transformation  of  the  estimated  optical  density.  To  do  this,  however,  requires  a knowledge 
of  the  conditional  variance  of  the  optical  density  in  the  two-dimensional  context,  which  is  accu- 
rate to  first  order  in  the  perturbations.  The  two-dimensional  image's  conditional  mean  has 
been  inferred  here  by  the  cascading  of  one-dimensional  conditional  means,  which  was  justified 
on  the  basis  of  the  relation  between  these  two  types  of  conditional  means  in  the  corresponding 
linear  case  without  perturbations.  Yet  it  is  difficult  to  see  how  this  reasoning  could  be  extended 
to  the  conditional  variance  perturbations,  because  they  form  part  of  the  output  of  the  nonlinear 
filter  but  not  of  the  corresponding  linear  one.  This  information  could  be  obtained,  however, 
through  a full  two-dimensional  formulation  of  the  problem,  such  as  that  described  in  Section 
II-C,  thus  enabling  the  conditional  mean  of  the  pixel  intensities  to  eventually  be  determined. 

Such  a formulation  would  be  much  more  formidable,  however,  and  is  beyond  the  scope  of  this 
report. 
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